n = size(A, 1);
f = rand(n, 1);
x = A \ f;
epsilon = 0.7;

D = diag(diag(A), 0);
R = A - D;
Q = R;
Q(abs(Q)<epsilon) = 0;

xs = zeros(n, 1);

for i = 1 : 10
    xs = D \ (f - Q * xs);
    norm(abs(x - xs))
end

% for i = 1 : 10
%     Q = R;
%     Q(abs(Q)<epsilon) = 0;
%     max(abs(eig(inv(D)*Q)))
%     epsilon = 2 * epsilon
%     figure
%     spy(Q)
% end
